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Two-component equal-mass Fermi gases, in which unlike atoms interact through a short-range 
two-body potential and like atoms do not interact, are stable even when the interspecies s-wave 
scattering length becomes infinitely large. Solving the many-body Schrodinger equation within a 
hyperspherical framework and by Monte Carlo techniques, this paper investigates how the prop- 
erties of trapped two-component gases change if a third or fourth component are added. If all 
interspecies scattering lengths are equal and negative, our calculations suggest that both three- and 
four-component Fermi gases become unstable for a certain critical set of parameters. The relevant 
length scale associated with the collapse is set by the interspecies scattering length and we argue that 
the collapse is, similar to the collapse of an attractive trapped Bose gas, a many-body phenomenon. 
Furthermore, we consider a three-component Fermi gas in which two interspecies scattering lengths 
are negative while the other interspecies scattering length is zero. In this case, the stability of the 
Fermi system is predicted to depend appreciably on the range of the underlying two-body poten- 
tial. We find parameter combinations for which the system appears to become unstable for a finite 
negative scattering length and parameter combinations for which the system appears to be made 
up of weakly-bound trimers that consist of one fermion of each species. 

PACS numbers: 



I. INTRODUCTION 

Over the past decade or so, the field of ultracold gases 
has seen tremendous breakthroughs. After reaching de- 
generacy in Bose [IU2|1 and Fermi (3| gases, the realization 
of an atom laser [JH) @] > of the Mott-insulator transi- 
tion 0] , and of the conversion from an atomic to a molec- 
ular gas [8[ followed. Many of the present-day studies 
take advantage of the tunability of the atom-atom scat- 
tering length in the vicinity of a so-called Fano-Feshbach 
resonance 

0,113 ■ As an external magnetic field is tuned 
through its resonance value, the sign of the scattering 
length changes [HI, Exactly on resonance, the scat- 
tering length is infinitely large, allowing for the study of 
strongly-correlated systems. Experimentally, the speed 
of the magnetic field ramp can be changed, allowing adi- 
abatic ramps, for example, and the ramp itself can be 
reversed. It is this versatility that made possible the 
experimental study of the BCS-BEC crossover, using ul- 
tracold atomic Fermi gases trapped in two different hy- 
perfine states [I3.[l4|. 

Using present-day technology, the realization of de- 
generate multi-component atomic Fermi gases appears 
possible in principle. The occupation of more than two 
different hyperfine states of the same species requires, 
neglecting for the moment possible losses, only moder- 
ate changes of current set-ups. A particularly promising 
candidate appears to be 6 Li and the coexistence of 
three hyperfine states has already been demonstrated for 
40 K [l6[ . Alternatively, a number of goups are presently 
pursuing the simultaneous trapping of three different 
atomic species [UIHI- I n the former scenario, the atomic 
masses of all components are equal, whereas in the latter 



scenario, the atomic masses of the different components 
differ. In either of these realizations of multi-component 
Fermi gases, all or some of the interspecies scattering 
lengths may be tunable thanks to the possible existence 
of magnetic or optical Fano-Feshbach resonances. This 
may open the possibility to experimentally investigate 
the stability and to study, provided an extended stable 
regime exists, the behaviors of multi-component Fermi 
gases as a function of the s-wave scatterig length. 

Using two different theoretical frameworks, this paper 
considers three- and four-component Fermi gases, and 
compares their behaviors with those of two-component 
Fermi gases. In particular, we ask how the stability 
of two-component Fermi gases changes when a third or 
fourth component are added. It is now well established 
that trapped two-component Fermi gases are stable even 
when the interspecies s-wave scattering length is negative 
and its magnitude is infinitely large [T^, [20, [U, [22|, [HI . 
The stability of inhomogeneous as well as of homogeneous 
two-component Fermi gases with attractive short-range 
interactions and arbitrary interspecies scattering length 
can be attributed to the Pauli exclusion principle (also 
referred to as Fermi pressure) , which introduces effective 
repulsive intraspecies interactions that more than com- 
pensate the attractive interspecies interactions. In con- 
trast, homogeneous Bose gases with negative scattering 
lengths are unstable; they can, however, be stabilized by 
an external confining potential as long as the product of 
the number of bosons and the s-wave scattering length is 
less negative than a certain critical value [U HH, Hy] . 

Section [TT] introduces the Hamiltonian of trapped 
multi-component Fermi gases as well as a simple "count- 
ing argument" that turns out to be quite useful in under- 
standing the stability of multi-component Fermi gases. 



2 



Section IIIII investigates the stability of three- and four- 
component equal-mass Fermi gases within a hyperspher- 
ical framework, focussing on the large and small parti- 
cle number limits. The physical picture emerging from 
the hyperspherical treatment is further investigated in 
Sec. IIV1 which solves the many-body Schrodinger equa- 
tion for short-range interactions using Monte Carlo tech- 
niques. Finally, Sec. [V] summarizes and connects our 
results with those available in the literature [2(J H3, [H, 
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II. HAMILTONIAN AND GENERAL 
CONSIDERATIONS 

The Hamiltonian H for an atomic Fermi gas with x 
components under external spherically symmetric har- 
monic confinement is given by 



x Na ( h 2 l 



a=l i=l 



2m a 

a</3 i=l j = l 



Here, the number of atoms of the ath component is de- 
noted by N a , and the total number of atoms is given by 
N, N = J2a=i Na- In Eq. ([1]), r ai denotes the position 
vector of the ith atom of the ath component, measured 
with respect to the center of the trap, and m a and uj a re- 
spectively the atomic mass and the angular frequency of 
the ath component. The potential V a p describes the in- 
teraction between an atom of the ath and an atom of the 
/3th component. This work considers a zero-range poten- 
tial (see towards the end of this section and Sec. and 
a purely attractive short-range potential (see Sec. lIVp . In 
both cases, we characterize the strengths of the V a p by 
the s-wave scattering lengths a a p. We assume that the 
two-body interactions are independent of spin, implying 
that the number of atoms in each spin state is conserved. 
Throughout this work, like atoms are taken to be non- 
interacting, implying a aa = for a — 1, • • • , X- 

In the most general case, the Hamiltonian given in 
Eq. (TTJ has x different m a , lu u and N a (a — 1, • • • , x), 
and x(x— 1)/2 different V a p (a, (5 = 1, • • • , x an d ot ^ (3), 
resulting in a tremendously large parameter space. To 
reduce the parameter space, we first consider the case 
where all m a , all o> Q , all N a , and all a Q( a (a / /3) are 
equal. A four-component gas of this type could, e.g., 
be realized by equally populating and trapping the four 
spin states of a fermionic atom whose ground state has 
vanishing total electronic angular momentum J but a 
non- vanishing nuclear spin / of 3/2. In this case, the 
scattering lengths between the different spin substates 
are equal and s-wave scattering between two atoms in 
the same spin substate are forbidden by symmetry. 

Before solving the Schrodinger equation for the Hamil- 
tonian given in Eq. |Tj), we present a simple counting 



TABLE I: Number N a tt of attractive interactions, number 
N rep of effectively repulsive interactions and ratio N rep /N a tt 
for finite and infinite N for a x-component Fermi gas (x = 2 
through 4) in which all interspecies interactions are equal (or 
resonant). 





X=2 X=3 X=4 


Natt 

N rev /N a tt (N finite) 
N rep /Natt (N -> oo) 


I TV 2 -A^ 2 -N 2 
f(f~l) f(f-l) f(f-l) 

JV-2 JV-3 JV-4 
JV 2JV 3JV 

1 1 1 

± 2 3 



analysis. The number N a tt of attractive pair interactions 
V a p, where again a ^ f3, is given by 



N a tt = 



N 2 x-l 
2 X 



C2) 



and the number N rep of effectively repulsive interactions, 
i.e., the number of like fermion pairs, by 

N N- 



X 



X 



(3) 



Table [J summarizes the values of N att , N rep and 
Nrep/N a tt for x = 2 , 3 and 4. The ratio N rep /N a tt (re- 
ported in the fourth and fifth row of Table U for finite 
and infinite TV, respectively) decreases with increasing 
X and approaches l/(x — 1) in the large N limit, indi- 
cating that the Fermi pressure becomes less important 
compared to the interspecies interactions as x increases. 
Another interesting scenario arises when each component 
is occupied by exactly one fermion, i.e., when x = N. 
In this case, no effectively repulsive interactions exist, 
i.e., N rep /N a tt = 0, and the system's ground state is the 
same as that of the corresponding TV-boson system. As 
pointed out already in the introduction, Bose gases with 
negative scattering length become unstable in the limit 
that the absolute value of the scattering length becomes 
large. This, together with the fact that two-component 
Fermi gases are stable for all scattering lengths, suggests 
that there exists a critical x~value beyond which multi- 
component Fermi gases with large negative interspecies 
scattering length are unstable. Sections IIIII and HVl show 
that Xcr = 3. 

In addition to Fermi gases in which all interspecies in- 
teractions V a f3 (a,/3 = 1, ■ • • , X an d ol ^ (3) are equal, 
we consider the scenario in which only a subset of inter- 
species interactions are "turned on" . In particular, we 
consider x-component Fermi gases with x ~ 1 equal and 
non-zero (or resonant) a a p. For the three-component 
system, we take an = 023 and 031 = 0. The number 
of attractive interactions is in this case by a factor of 
2/3 smaller than in the case with three resonant inter- 
actions, thus increasing the ratio of N rep /N a tt from 1/2 
to 3/4 in the large TV limit. For the four-component sys- 
tem, two different "non-trivial" possibilities for turning 
on only x ~ 1 interactions exist: (i) 012 = 013 = an 
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and 023 = 034 = 0,24 = 0, and (ii) a\2 = 023 = 034 and 
a i3 = a i4 = a 24 = [the configuration an — 023 = 034 
and (Z13 = &i4 = a.24 = 0, e.g., is equivalent to (i); the 
configuration 012 = 023 = 031 and 014 = 024 = 034 = 0, 
e.g., is trivial in the sense that it can be broken up into 
a three-component system with all resonant interactions 
and a single non- interacting component]. For the non- 
trivial configurations, the number of attractive interac- 
tions is half as large for the case of three resonant inter- 
actions as in the case of six resonant interactions, thus 
increasing the ratio of N rep /N att from 1/3 to 2/3 in the 
large N limit. This counting analysis indicates that the 
values of the ratio N rep /N att for large N for three- and 
four-component gases with x ~ 1 resonant a a p are be- 
tween those for two- and three-component Fermi gases 
with x(x ~ resonant a a p. Thus, the question arises 
whether multi-component Fermi gases with \ ~ 1 reso- 
nant interactions are stable for all scattering lengths and 
N values, or whether they become unstable for a certain 
critical parameter combination. 

To analyze the stability of Fermi gases quantitatively, 
one may attempt to describe the interspecies atom- 
atom interactions by a zero-range Fermi pseudopotential 

v s (f) m, 

Tr 2-Kh 2 a s , A . 

Vsif) = S(r), (4) 

A 4 

which is directly proportional to the s-wave scattering 
length a s . Here, fx denotes the reduced mass of the two 
interacting atoms. This pseudopotential has been em- 
ployed successfully to predict many properties of dilute 
Bose gases. For example, the interaction potential given 
in Eq. ^ together with a Hartree wave function cor- 
rectly predicts that tra ppe d Bose gases with negative a s 
become unstable if [24l |42|| 

(7V_i)^<_ .575, (5) 

a-ho 

where a; lo denotes the oscillator length, a/, = 
yjh/ {2^uj). In general, the true atom-atom interaction 
can be replaced in the long- wavelength limit by the pseu- 
dopotential given in Eq. (|4]) provided the system is di- 
lute, i.e., if n(0)|a s | 3 <§C 1, where n(0) denotes the peak 
density. Assuming N is not too small, one finds that 
n(0)|a s | 3 is much smaller than one when (N— l)a s equals 
— 0.575a^ o ; consequently, the pseudopotential predicts 
the instability of dilute Bose gases correctly. 

Applied to two-component Fermi gases, the bare pseu- 
dopotential employed within a hyperspherical framework 
predicts that the system becomes unstable if [43| 

M°K < -1.22, (6) 

where kp(0) denotes the noninteracting Fermi wave vec- 
tor at the trap center. (A mean-field analysis predicts 
a slightly more negative critical value of kp(0)a s < 
— 7r/2 [3, |45|, S^|.) However, at the critical parameter 
combination fci?(0)a s = —1.22, the diluteness criterium, 



which can be written as (&)F(0)|a s |) 3 <C 1 for fermions, 
is violated and the instability prediction, Eq. ([S]), can- 
not be trusted. Indeed, Eq. © disagrees with the 
experimental findin g th at two-component Fermi gases 
are stable [3 Q M S ES Efll- To resolve this 
disagreement, two independent studies recently intro- 
duced density-dependent "renormalization schemes" of 
the scattering length a s entering into Eq. (UJ) [5l|, 
Using these "renormalization schemes" , the (modified) 
density-dependent pseudopotential can be applied to de- 
scribe strongly-interacting Fermi gases and predicts, in 
agreement with experimental [l3l. fl4l |47l |H, [4^, [5(| and 
other theoretical [ijj [2(| I2H l22l loll] works, that two- 
component Fermi gases are stable even in the strongly- 
interacting unitary regime. The fact that the Fermi pseu- 
dopotential has to be modified when applied to fermions 
was already suggested earlier (see, e.g., Refs. 0, 
and is well known in the nuclear physics community (see, 
e.g., Refs. [55l 1561] and references therein). 

The instability prediction given in Eq. © for two- 
component Fermi gases can be readily generalized to 
multi-component Fermi gases with x(x — l)/2 and with 
X — 1 resonant interactions (in the case of % — 1 resonant 
interactions, we consider the non-trivial scenario intro- 
duced above, in which none of the components can be 
separated off). The right hand side of Eq. ^ has to 
be multiplied by l/(x — 1) m the former case, and by 
x/[2(x — 1)] in the latter case. Since the diluteness cri- 
terium is fullfilled approximately in the all resonant case 
at the predicted collapse point [e.g., (kp(0)a s ) 3 w —0.23 
and —0.067 for x = 3 and 4, respectively], the prediction 
that multi-component Fermi gases with all resonant in- 
teractions become unstable for a finite |a s |, derived using 
the "bare" interaction given in Eq. ((4]), may in fact be 
correct. Sections IIIII and IIVI confirm this. In the case 
with x ~ 1 resonant interactions, the bare Fermi pseu- 
dopotential, Eq. , predicts that the collapse occurs at 
(A:F(0)a s ) 3 w —0.76 and —0.54 for x = 3 and 4, respec- 
tively. This suggests that the bare pseudopotential can- 
not be used and that the instability prediction derived 
using it may not be correct (see Sec. IIVI for a many-body 
analysis). 



III. HYPERSPHERICAL FRAMEWORK 

A. iV-particle system 

This section investigates the stability of three- and 
four-component Fermi gases within a hyperspherical 
framework [3?], IH, H3, HH . Throughout this section, we 
assume that all angular trapping frequencies are equal, 
i.e., uj a — uj for a = 1, • • • , x- To gain insight into the sys- 
tem's behavior, we employ hyperspherical coordinates: 
the 3iV coordinates are divided into a hyperradius R and 
3N — 1 hyperangles, collectively denoted by f2 [5^ ]. In 
the following, only the definition of the hyperradius R, 
which can be thought of as a measure of the size of the 
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system, is needed, 



X N a 



R = U/ EE m «^ 



(7) 



OL — l i—1 



As before, the position vectors r a i are measured with re- 
spect to the center of the trap. In Eq. 0, M denotes 

the total mass of the system, i.e., M = J2a=i 2i=i m a- 
Using these coordinates, the many-body wave function 
^(R, 0) can be expanded in terms of a set of ^-dependent 
channel functions & u (fl;R), which depend parametri- 
cally on R, and i?-dependent weight functions F un (R), 



(8) 



9{R, n) = J2 R {1 ' 3N)/2 F m (R)<P,(Q; R). 



In the adiabatic hyperspherical approximation 57, 60, 
l6ll . l62l . f&3 | , which neglects off-diagonal coupling elements 
and additionally restricts the sum in Eq. ([5]) to a single 
term, the solution of the many-body Schrodinger equa- 
tion reduces to determining an i?-dependent effective po- 
tential V V (R), which includes part of the kinetic energy as 
well as the effects of the two-body interactions, and then 
solving an effective one-dimensional radial Schrodinger 
equation in the hyperradius R, 



2M dR 2 



+ Vu(R) + V tra p(R) 



F un {R) = 



where 



V trap {R) = -Muj 2 R 2 . 



(9) 



(10) 



In general, the calculation of the effective potentials 
V V {R) is, at least for many-body systems, just as hard 
as solving the many-body Schrodinger equation itself. In 
certain circumstances, however, the effective potentials 
or their functional form can be obtained more easily. 

For weakly-interacting dilute equal-mass two- 
component Fermi gases, e.g., the effective hyperradial 
potential Vq(R) consists of two parts [13, EH]: The first 
part, which is proportional to Ckin/R 2 , arises from the 
kinetic energy operator, and the second part, which 
is proportional to Ci nt /R 3 , accounts for the particle- 
particle interactions. Ckin is positive, and Ci n t is directly 
proportional to the s-wave scattering length a s |43[. If 
\a s \ is sufficiently small (<z s < 0), the small R region, 
where the attractive Ci n t/R 3 term dominates over the 
repulsive Ckm/R 2 term, is separated by a "barrier" 
from the R region where the repulsive Ckin/R 2 term 
dominates over the attractive Ci„t / R? term. This barrier 
prevents the Fermi gas from collapsing to a cluster-like 
many-body bound state, thus explaining the stability 
of weakly-interacting dilute equal-mass two-component 
Fermi gases (see also Ref. (5l|). 

Analytical expressions for V V (R) for equal-mass two- 
component Fermi gases have also been obtained in the 



strongly- interacting limit [37|, [6J|. For infinitely large 
interspecies scattering lengths a a p (a ^ (3), the adia- 
batic approximation introduced above is — for a class of 
states that arise assuming boundary conditions consis- 
tent with contact interactions when the distance between 
each pair of particles goes to zero [64( , referred to as uni- 
versal states in the following — exact and the universal 
effective potentials V V {R) are given by 



V V {R) = 



h 2 (i + c v ) 

2MR 2 



(11) 



In Eq. (Tnj) , the C u denote R- independent constants that 
arise from the integration over the 3iV — 1 hyperangles. 
The functional form given in Eq. (fTTj) has also been de- 
rived by explicitly solving the many-body Schrodinger 
equation for a class of variational wave functions using 
hyperspherical coordinates [37| . Using dimensional argu- 
ments, Eq. (jTTJ) can be understood as follows: The term 
of V V (R) that accounts for the particle-particle interac- 
tions has — because of the absence of any other length 
scale in the problem — to scale in the same way with R 
as the 1/R 2 term that accounts for the kinetic energy 
contribution. The total universal effective hyperradial 
potential curves at unitarity thus have a simple func- 
tional form: V V (R) dominates at small R and approaches 
plus or minus infinity in the R — > limit, depending on 
whether the quantity 1 + C v is positive or negative, and 
Vtrap(R) dominates at large R. 

For a two-component unitary Fermi gas, a number of 
C values are known. For N < 6, selected C„ with v < 2 
have been obtained using a correlated Gaussian (CG) 
basis set expansion- type approach [65[. For N < 30, 
upper bounds for the Co have been obtained by the FN- 
DMC method 

In the large N limit, the value of Co for equal-mass 
two-component Fermi gases at unitarity can be related 
to the universal parameter j3 of the homogeneous sys- 
tem using the hyperspherical framework of Ref . 13711 that 
employs a density-dependent scattering length [521 ] , lead- 
ing to Co = (3 [37]. Alternatively, this relationship can 
be derived by applying the local density approximation 
to the trapped system (see, e.g., Ref. [65J). It is gener- 
ally believed that the most accurate value for (3 has been 
obtained by the FN-DMC method [Hill, = -0.58. 
Since l+f3 > 0, the universal hyperradial potential curve 
Vq(R) + Vt ra p(R) for large N (shown by a solid line in 
Fig- [Busing (3 — —0.58) is repulsive at small R, prevent- 
ing the collapse of the two-component unitary Fermi gas 
towards cluster-like bound states. The hyperspherical 
potential curve picture reveals an intuitive way of un- 
derstanding the stability of the energetically lowest-lying 
gas-like state of two-component unitary Fermi gases. 

In Fig. [TJ the hyperradial potential V(R) and the hy- 
perradius R are scaled by E^j and Rni, respectively, 
where Eni denotes the energy of the non-interacting 
^-component Fermi system and Rni the square root 
of the expectation value of R 2 calculated for the non- 
interacting ^-component Fermi system, i.e., Rni = 
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FIG. 1: (Color online) Hyperradial potential curve Vo(R) + 
Vtrap(R) as a function of the hyperradius R for \ = 2 (solid 
line) , x — 3 (dotted line) and \ — 4 (dashed line) in the large 
N limit. All interactions between unlike atoms are character- 
ized by an infinite scattering length, and the coefficient Co is 
taken to be (x~ with /3 = —0.58. Both length and energy 
are scaled by the corresponding values of the non-interacting 
system (see text). 



^(R 2 ) ni [33, Sa| ■ Using the virial theorem [64,l67| ; Rni 
can be related to Eni, 



Rni 



h 

~M~uj 



Em 



(12) 



In Fig. [IJ En i and Rni have been evaluated by writing 
En i as (X + 3N/2)huj and using the leading order of A for 
large closed-shell systems, A = (3iV) 4 / 3 /4. We note that 
the scaled quantity R/ Rni remains finite in the limit that 
N goes to infinity, allowing the thermodynamic limit to 
be taken. 

For all-resonant three- and four-component gases at 
unitarity, the functional form give n in Eq. (|II[) remains 
valid for universal states (37L l64l| . For these systems, 
the values of C v are not known (the only exception be- 
ing the N = 3 system [68|, |69|, see below). It has been 
argued [37| . parametrizing V a a by a density-dependent 
zero-range two-body potential [52j and neglecting the ex- 
istence of weakly-bound trimer or tetramer states, that 
the coefficient C is in the large N limit determined by 
the parameter [3 of the unitary two-component Fermi gas, 
i.e., Co = (x~ Due to the lack of benchmark results 
for all-resonant three- and four-component Fermi gases, 
there is some ambiguity in how the simplest adaptation 
of the renormalization scheme, originally designed to de- 
scribe the physics of two-component Fermi gases [52j . 
is applied to three- and four-component Fermi gases, 
thus introducing some uncertainty about the exact values 
of Co that describe unitary three- and four-component 
Fermi gases. Using C = (%-l)/3 and /3 = -0.58 [lHH, 
(I + Co) is negative for three- and four-component gases, 
(1 + Co) = -0.16 and (I + C ) = -0.74 for x = 3 and 



4, respectively [TO] ■ The resulting hyperradial potential 
curves Vo(R) + V tra p{R) are shown in Fig. Q] in the large 
N limit by dotted and dashed lines for x = 3 and 4, 
respectively. The attractive small- R behavior is due to 
the negative values of (1 + C ), and suggests that unitary 
three- and four-component Fermi gases with all resonant 
interations are unstable. Note that the coefficient Co and 
N a tt/N rep both scale as 1 — x with x, reflecting the fact 
that Vo(R) is in part determined by the pair interactions. 

Interestingly, both the bare pseudo-potential (see the 
discussion towards the end of Sec. [TTJ) and the density- 
dependent pseudo-potential (used in this section within 
the hyperspherical framework) predict that x-component 
Fermi gases, \ > 2, with all resonant interactions become 
unstable for a finite value of the interspecies scattering 
length. The critical value predicted by the bare inter- 
action is presumably not negative enough while that for 
the density-dependent interaction would presumably be 
too negative (we expect that the renormalization scheme 
originally developed for the two-component Fermi gas 
"over-renormalizes" the scattering length). 

The analysis above can be extended to the case where 
only x — 1 interspecies scattering lengths are infinite; as 
before, we focuss on what we defined in Sec. [IT] as "non- 
trivial" scenarios. In the case of x — 1 resonant inter- 
actions, the bare pseudopotential is expected to fail (see 
Sec. [TTJ) . Using, just as above, density-dependent zero- 
range two-body interactions to describe the V a p with 
non-zero a a p, the coefficient Co becomes in the large N 
limit ^(x~ just as in the all-resonant case, Co scales 
with x in the same way as N at t/N rep . Using [3 — —0.58, 
(I +Co) is positive for x — 3 and 4 (the corresponding po- 
tential curves are shown in Fig. [21) ; the repulsive small- R 
behavior suggests that three- and four-component uni- 
tary Fermi gases with x — 1 resonant interactions are 
stable. We note that the analysis in hyperspherical coor- 
dinates, which employs density-dependent interactions, 
has two short-comings: i) The derivation of the effec- 
tive potential curves assumes the existence of a single 
length scale, the hyperradius R; since some of the inter- 
species scattering length are zero while others are res- 
onant, a more accurate description would presumably 
introduce an additional length scale and be based on 
a "hyperradial vector" as opposed to a single hyperra- 
dius. ii) The derivation neglects the possible existence of 
weakly-bound trimer (see the next section) or tetramer 
states, or said differently, the existence of non-universal 
states. Section IIVBI shows that the behaviors of multi- 
component Fermi gases with x > 2 may differ depending 
on whether or not such states exist. 



B. Three-particle system 

This section discusses the behaviors of three distin- 
guishable fermions interacting through zero-range poten- 
tials with equal masses and equal trapping frequencies. 
For these systems, the effective hyperradial potential 
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FIG. 2: (Color online) Hyperradial potential curve Vo(R) + 
Vtrap(R) as a function of the hyperradius R for \ — 3 (dotted 
line) and x = 4 (dashed line) in the large N limit for \ — 1 res- 
onant interactions (the scattering lengths a\p, j3 = 2, • • • 
are infinite and all other scattering lengths are zero). The 
coefficient Co is taken to be — (x — with (3 — —0.58. Both 
length and energy are scaled* by the corresponding values of 
the non- interacting system (see text). The plotting range of 
the axis is the same as in Fig. [T] for ease of comparison. 
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FIG. 3: (Color online) Solid, dash-dotted and dashed lines 
show the hyperradial potential curve Vo(R) + Vtrap(R) as a 
function of the hyperradius R, scaled by Rni, for three dis- 
tinguishable fermions interacting through zero-range interac- 
tions with a 3 — — O.ldho, — 0.15a/, o and — Q.2aho, respectively. 
Note that R is defined without the center-of-mass motion, 
implying Emi = 3hu>. 



curves V V {R) can be obtained analytically for any com- 
bination of scattering lengths an , 023 and 0,31 [71], uM ■ 

Figure [3] shows Vq(R) + V tra p{R) for the three-particle 
system as a function of R for three different negative 
scattering lengths a s , i.e., a s = — 0.1a/ lo , — 0.15a^ o and 
— 0.2a,ho (here, a s — a a p with a ^ (3) as a function of 
the hyperradius R. Throughout this section, the hyper- 



radius R is defined without the center-of-mass motion 
and scaled by the hyperradius Rni of the non-interacting 
system without the center-of-mass motion. For the dis- 
cussion that follows, the difference between the hyper- 
radius defined without the center-of-mass and that de- 
fined in Eq. (J7J) is irrelevant. For the scattering lengths 
shown, the hyperradial potentials show a barrier around 
R w 0.2Rni (i.e., R w 0.3a; lo ), which separates the large 
R region where gas-like states live from the small R re- 
gion. As a s becomes more negative, the barrier height 
decreases and moves to slightly larger R values. The 
critical scattering length a c at which the barrier disap- 
pears (a c « —0.39aho) provides a first estimate for when 
the three-body system becomes "unstable" against the 
formation of tightly-bound trimer states with negative 
energy, assuming such states exist; the conditions for the 
existence of negative energy states are discussed below. 
A more accurate estimate would account for the zero- 
point energy of the quantum system, resulting in a some- 
what less negative critical scattering length that is in fair 
agreement with the GP prediction given in Eq. © (see 
Ref. for a discussion of the iV-boson system). The 
lifetime of the gas-like three-body state can be estimated 
by calculating the tunneling rate through the potential 
barrier as a function of a s [571 ] . When the barrier height is 
large compared to the energy of the non-interacting sys- 
tem (small \a s \), the lifetime of the gas-like trimer state 
is much larger than l/u>; in this case, the system can be 
considered stable (it is, in fact, metastable). The tun- 
neling rate is appreciable only when the barrier height 
becomes comparable to the energy of the non-interacting 
system, i.e., when a s approaches a c . 

To obtain the energy spectrum in the adiabatic approx- 
imation, we solve the one-dimensional radial Schrodinger 
equation [Eq. ©]. To this end, we add a term Vsr(R) 
that is repulsive for R < R c and negligibly small for R > 
R c to V V (R), Vsr(R) = 2(1 - V3tanh 2 (i?/i? c ))/(2mi? 2 ). 
Vsr(R) cures the divergence of the V v (R) at small R, 
which is related to the Thomas collapse [73]. When the 
barrier height is large compared to the energy of the non- 
interacting system and the short-range length scale R c 
is chosen sufficiently small, this "ad-hoc" regularization 
should affect the eigenenergies of the gas-like states at 
most weakly. The eigenenergies of molecular-like states 
(states living in the small R region), however, should de- 
pend comparatively strongly on the particular value of 
R c . Thus, some properties of the system are expected 
to be non-universal. Indeed, our calculations for differ- 
ent R c show that the first state with negative energy 
appears when the scattering length a s is approximately 
equal to — 8AR C . For R c = 0.001a/i O , e.g., the first neg- 
ative energy state appears for a s — 0.0084a^ o , a value 
much less negative than a c . At this scattering length, the 
barrier height is large compared to the energy of the non- 
interacting system and the gas-like system is — because 
of the existence of a negative trimer state — metastable. 
If R c is chosen so that the absolute value of the scattering 
length at which the first bound state appears is compa- 
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rable to or larger than \a c \, then the gas-like state would 
be the true ground state of the system for \a s \ < \a c \. For 
realistic dilute alkali gases with all-resonant interactions, 
we expect that three-body bound states appear for scat- 
tering lengths less negative than a c ; consequently, trimcr 
as well as larger three-component systems with a s > a c 
are expected to be metastable. The picture developed 
here within the adiabatic approximation remains qualita- 
tively valid if the coupling between channels is included. 

In addition to the all-resonant three-particle system, 
we consider the trapped three-particle system with two 
resonant interactions (a s = ayi = 023 and 031 = 0). In 
this case, the hyperradial potential barrier disappears at 
a more negative scattering length than in the all-resonant 
case, i.e., at a c » — 1.00<z/j O . An analysis of the tunneling 
probability suggests that the lifetime of the two-resonant 
system is, considering the same <z s , enhanced by a factor 
of about ten compared to the all-resonant system. A life- 
time enhancement is expected since the ratio N at t/N rep 
(see Sec. |TTJ) is smaller for the three-component system 
with two resonant interactions than for that with three 
resonant interactions. 

Calculating the trimer energies in the adiabatic ap- 
proximation, we find that the first three-body state with 
negative energy appears at a s sa — 2500i? c . For exam- 
ple, a negative energy state appears for a s ps — 2.5a/j if 
R c = O.OOla^ and for a s « -0.25a ho if R c = O.OOOlaw 
In the latter case, the disappearing of the hyperradial 
potential barrier is accompanied by a "collapse" of the 
metastable gas-like state to a cluster-like state. In the 
former case, in contrast, the disappearing of the hyperra- 
dial potential barrier is accompanied by the lowest-lying 
gas-like state evolving smoothly to a cluster-like state 
with negative energy. Although the short-range param- 
eter R c cannot be straightforwardly related to the pa- 
rameters of typical atom-atom potentials, it seems plau- 
sible that realistic three-component Fermi systems with 
two resonant interactions could fall into either of the two 
categories discussed here. Section HVl investigates the be- 
haviors of three-component many-body systems and at- 
tempts to determine how the system's behaviors depend 
on the range of the underlying finite-range two-body po- 
tential (this range can, roughly speaking, be connected 
with the short-range parameter R c ). 

We note that the negative scattering length a s at 
which the first weakly-bound trimer state appears (as- 
suming no deeply-bound states exist Vcan be estimated 
using the following "rule of thumb" [3, [lE [ill : |a 8 | = 
R c exp(7r/so)i where sq is the coefficient that determines 
the energy spectrum of the lowest hyperradial potential 
curve at unitarity (so = 1.006 for a s = an = 023 = 031 , 
and So = 0.413 for a s = <zi2 = &23 and 031 = [68j| ) . 
This rule of thumb predicts a negative scattering length 
at which the first weakly-bound trimer state appears of 
— 23i? c and — 2000i? c for three and two resonant interac- 
tions, respectively. These values are to be compared with 
— 8.4i? c and — 2500i? c found numerically (see above). 



IV. MONTE CARLO TREATMENT 

This section discusses the solutions of the many-body 
Schrodinger equation for the Hamiltonian given in Eq. ([1]) 
obtained by Monte Carlo techniques as a function of the 
scattering length. Section [IV Al introduces the VMC and 
FN-DMC methods, and Sec. IIV Bl presents our results. 



A. Variational and fixed-node diffusion Monte 
Carlo method 

Using the VMC and FN-DMC methods [h| , we deter- 
mine solutions of the many-body Schrodinger equation 
for two different purely attractive short-range model po- 
tentials V a p: a square well interaction potential r , 



V a T(r) = 



-V a i3,o for r < R a/ 3 : o 
for r > R a/ 3fi, 



(13) 



and a Gaussian interaction potential 



0' 



-V a p fi exp 



2R2 a0,O J 



(14) 



Both potentials depend on a depth V a p,o, V a p,Q > 0, 
and a range R a p.Q. In our calculations, we set all R a f3,o 
to the same value, R a pfi <C a^ , and adjust the depths 
Vq/3,0 until the s-wave scattering lengths a a p assume the 
desired values. In Sec. IIVB1 V a pfi and R a p,o ar e chosen 
so that the potential V a p supports no two-body s-wave 
bound state and so that a a p < 0. The use of two different 
functional forms for the interaction potentials V a p allows 
for exploring the dependence of the results on the details 
of the short-range physics. 

We now discuss how the solutions of the many-body 
Schrodinger equation are obtained by the variational 
Monte Carlo method. The functional form of the vari- 
ational wave function ipx is chosen on physical grounds 
(see below), and ipr is parametrized in terms of a set of 
parameters p, which are optimized so as to minimize the 
energy expectation value E{p) (77| . 



(4> t (p)\h\Mp)) 
(ip T (p)\ip T (p)) 



(15) 



Since i/'t depends in general on the 3N coordinates of the 
system, the integration on the right hand side of Eq. (fT5"|) 
is high-dimensional; this implies that standard tech- 
niques such as the Simpson rule [78| |. which scale, roughly 
speaking, exponentially with the number of degrees of 
freedom cannot be used for its evaluation. Instead, we 
evaluate the high-dimensional integral by Monte Carlo 
techniques using standard Metropolis sampling [7^] . The 
stochastic evaluation of the integral introduces a statis- 
tical uncertainty, which can be reduced by increasing the 
number of samples K included in the evaluation of the 
expectation value (the errorbar decreases as 1/y/K with 
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increasing K). We denote the energy expectation value 
obtained by the VMC method by E VMC . 

The variational wave function ipx for multi-component 
Fermi gases is written as (see Refs. [2l|, [22|, HH, El, H3, 
HHH^I for Monte Carlo studies of two-component systems 
and Ref. [35[ for a Monte Carlo study of three-component 
systems) 

^ t = n n n f^Ki - x n^ Q (r Qi ) 

a<f3 i— 1 j — 1 a—li—1 

X 

x Y[ A{(f>l (f al ), ■ ■ ■ ,<t>N a (r a N a )), (16) 

a=l 

where A denotes the anti-symmetrizer and f a p, (f a and 
4>i denote two- or one-body functions whose functional 
form is discussed in detail in the following. 

The f a p denote two-body correlation factors. For 
small r, f a p coincides with the zero-energy scatter- 
ing wave function for two particles interacting through 
V a p. Beyond some matching value R a /3, m , which is 
treated as a variational parameter, f a p is taken to be 
Cq/3,i + c a/ 3 t 2 exp(— c a /3r), where c a p denotes a variational 
parameter and c a p_i and c a /3^ are determined by the 
condition that f a p and its derivative are continuous at 
r = R a f3, m - The R a /3, m and c a p are optimized under the 
constraint that f a p > for all r. 

The ip a denote Gaussian orbitals with widths b a , 
Pair) = ex P(~ r2 /(2^a))- The parameters b a control the 
size of the Fermi cloud, and are optimized variationally. 
In the non-interacting case, b a = a; lo . For a a p < 0, 
the system becomes more compact due to the attractive 
interactions, resulting in a smaller optimal value of b a . 
If all m a and all b a are the same, i.e., m a = m and 
b a = b, then the product n«=i fltTi fai^ai) reduces to 
exp(—NR 2 / (2b 2 )) and the variational parameter b is pro- 
portional to the mean hyperradius of the variational wave 
function. 

The <f)i(r) are given by H nx (x)H ny (y)H nz (z), where 
(n x ,n y ,n z ) = (0,0,0) for i = 1, (n x ,n y ,n z ) = (1,0,0) 
for i — 2, (n x , n y , n z ) = (0, 1, 0) for i = 3 and so on, and 
the H n denote Hermite polynomials of degree n. 

The anti-symmetry and nodal surface of tpT are de- 
termined by the x Slater determinants [second line of 
Eq. |H])], which contain the one-body functions <fii. The 
nodal surface of -ipT coincides with that of the non- 
interacting multi-component Fermi gas. For small |a a ^|, 
o-a/3 < 0, this is expected to be a very good approxima- 
tion. It has been shown that the nodal surface used here 
also provides reasonably tight bounds for the energies of 
trapped two-component gases [H, HH, [H[ all the way to 
unitarity. Whether this holds true for multi-component 
Fermi gases with \ > 2 should be addressed in more de- 
tail in follow-up work. Note that the guiding function ipx 
used throughout this work does not build in any "pair- 
ing physics" from the outset (see Sec. [V] for further dis- 
cussion). The product <fi a (r)ipi(r), i — 1, ,N a , with 
b a = a,ho coincides with the ith (non-normalized) single- 



particle harmonic oscillator wave function. We choose 
to write t/'t in terms of the product fyaSPi instead of the 
harmonic oscillator functions themselves, because this al- 
lows the widths b a of the Gaussians to be varied without 
changing the nodal surface of ipT- 

In addition to the VMC method, we apply the FN- 
DMC method [z3,[13. The FN-DMC method uses the 
variational wave function tpT as a so-called guiding func- 
tion and determines the energy of a state whose nodal 
surface is identical to that of tpr- It can be shown that 
the FN-DMC energy provides an upper bound to the 
lowest eigenenergy of the eigenstate that has the same 
symmetry as tfjr Hl|. If the nodal surface of ipr coin- 
cides with that of the true eigenfunction, then the FN- 
DMC method results, within the statistical uncertainty, 
in the exact eigenenergy. In general, the quality of the 
FN-DMC energy depends crucially on the construction 
of the nodal surface of ipT- As in the VMC method, 
expectation values calculated by the FN-DMC method 
have a statistical uncertainty, which can be reduced by 
increasing the computational efforts. 



B. Monte Carlo results 

We first consider three- and four-component Fermi 
gases with one atom per component and equal masses 
and trapping frequencies. If all interspecies scattering 
lengths are equal, these Fermi gases are described by the 
same Hamiltonian as the corresponding Bose gas and the 
system should become unstable towards the formation of 
negative energy states, characterized by small interpar- 
ticle distances, when the inequality given in Eq. (0 is 
fulfilled. For N = 3 and 4, this implies critical scattering 
lengths of a s s» — 0.29a/j O and — 0.19a/ lo , respectively. 

To investigate how this instability or collapse arises 
within the many-body Monte Carlo framework, we per- 
form DMC calculations for the three- and four-particle 
systems interacting through a square well potential with 
Rq = O.Ola/,,0 (the subscripts a and f3 of R a /3fi have been 
dropped for notational convenience). The depth is cho- 
sen so that a s takes the desired value and so that the 
potential supports no two-body bound states. For small 
\a s \, the DMC method results in energies that agree very 
well with the mean-field Gross-Pitaevskii (GP) energy 
for both N = 3 and 4. As a s becomes more negative 
(of the order of — 0.06a,h, o [HI), the DMC walkers sample 
at first exclusively positive energy configurations; after a 
certain propagation time some walkers sample configu- 
rations with negative energy that correspond to tightly- 
bound states. The existence of tightly-bound three-body 
states, which depend on the details of the underlying two- 
body potential, for scattering lengths a s that are less neg- 
ative than the critical scattering lengths determined at 
the mean-field level or within the hyperspherical frame- 
work was already pointed out in Sec. IIII Bl Our calcula- 
tions here show that the DMC algorithm "sees" these 
three-body bound states (as well as four-body bound 
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FIG. 4: (Color online) Circles, squares and triangles show 
the variational energy Evmc for TV = 3 atoms interacting 
through a square well potential (range Ro = 0.01ajj O ) with 
a s = — Q.05a,ho, —Q.ldho and — 0.15a/i O , respectively, as a 
function of the variational parameter b [see the discussion 
following Eq. (|16p ], Evmc is scaled by the energy i?jv/ of the 
non-interacting system, Eni = 4.5ftc<J. Dotted lines connect 
data points for ease of viewing. 



states) if the simulation time is sufficiently long and \a s \ 
is sufficiently large. If |o s | is not too large, the DMC al- 
gorithm does — if the initial walker configurations corre- 
spond to gas-like states — not "know" about the tightly- 
bound states with negative energy. 

To gain further insight, we treat the N = 3 system by 
the VMC method. For a given a s , we fix the variational 
parameters R m and c (their exact values are not impor- 
tant for the discussion that follows) and vary b (we drop 
the subscripts a and /3 of R a /3, m and c a p, and the sub- 
script a of b a ). Circles, squares and triangles in Fig. [4] 
show the variational energy Evmc for N = 3 as a func- 
tion of b for a s — —0.05ah o , — O.ld/jo and — 0.15a/j O , re- 
spectively. For a s — — 0.05a/j O , Evmc increases if the 
Gaussian width b is larger than about aho (the system 
expands compared to its optimal size, thereby reducing 
the attraction felt between the atoms) and if b is smaller 
than about aho (the system shrinks compared to its op- 
timal size, thereby increasing the kinetic energy). For 
a s = — O.la/jo and — 0.15a/ lo (squares and triangles in 
Fig. ID), the variational energy assumes a local minimum 
at b w aho, shows a local maximum at b w 0.25a/j O for 
a s = — O.ldfco and at b w 0.5oft o for a s — — 0.15a/j O , and 
becomes negative for yet smaller b. We refer to the local 
maximum as an "energy barrier" that separates the local 
minimum at b « aho from a global minimum that exists 
at smaller b values. Finally, for more negative scattering 
lengths (not shown in Fig. [3]) no energy barrier exists for 
the variational wave functions considered. 

The energy barrier discussed here for small s-wave 
interacting Bose gases also exists in three-dimensional 
dipolar Bose gases [H[ and one-dimensional Bose 



gases [86| (the second part of the Appendix in Ref. |85j 
provides a detailed discussion). As in those earlier stud- 
ies, we interpret the existence of the energy barrier as 
an indication that the Hilbert space is divided into two 
nearly orthogonal spaces if a s is quite a bit less nega- 
tive than the critical scattering length predicted at the 
mean-field level: Gas-like states live in one region of the 
Hilbert space, and cluster-like bound states in the other. 
While this reasoning leads to an intuitive understand- 
ing of the stability and decay of Bose gases with neg- 
ative scattering length, the question arises why the en- 
ergy barrier disappears for scattering lengths a s that are 
notably less negative than the critical scattering length 
predicted by the GP equation. We believe that this is 
due to the existence of tightly-bound states with nega- 
tive energy. If \a s \ is sufficiently large, the overlap of the 
variational wave function with eigen functions of cluster- 
like bound states increases with decreasing b. For non- 
vanishing overlap, Evmc provides a rather poor upper 
bound to the true ground state energy of the system and 
not an upper bound to the energetically lowest-lying gas- 
like state. This implies that our VMC calculations do 
not allow for a reliable determination of the critical scat- 
tering length. Instead, they indicate that the existence 
of the energy barrier gives rise to the stability of the 
gas-like state and that the disappearing of the energy 
barrier qualitatively explains how the instability of Bose 
gases and multi-component Fermi gases with one atom 
per component arises when a s becomes too negative. 

We note that the picture developed here based on our 
VMC and DMC calculations is closely related to the anal- 
ysis based on the hyperspherical formulation presented in 
Sec. Hn] and in Refs. [33, M, M, H ■ In the VMC cal- 
culations, the role of the hyperradius R is played by the 
Gaussian width b [see discussion following Eq. ([Tfi|l ]. Fur- 
thermore, variational mean-field calculations have been 
interpreted in much the same way [89l . |9C| . 

We now investigate the behaviors of multi-component 
Fermi gases with more than one atom per component for 
which all interspecies scattering lengths a a p with a ^ (3 
are resonant; as before, we set a a p = a s . We first solve 
the many-body Schrodinger equation for a square well 
potential with range R a p,a = O.Olaho by the FN-DMC 
method. Circles in Fig. O show the FN-DMC energy for 
a four-component Fermi gas with four atoms per com- 
ponent (this corresponds to a closed shell) as a func- 
tion of the s-wave scattering length a s . The energy de- 
creases approximately linearly with decreasing a s , sug- 
gesting that this weakly-interacting system can be de- 
scribed to a good approximation perturbatively. Assum- 
ing zero-range interactions with scattering lengths a s and 
treating the system within first order perturbation the- 
ory, we find E « 36huj + ca s /aho, where c = 93hw / \/2tt , 
for N = 16. This perturbatively calculated energy (solid 
line in Fig. [5]) describes the FN-DMC energies (circles) 
very well. By additionally treating the N = 17 fermion 
system with Ni = 5 and N 2 = N 3 = N 4 = 4 (squares 
in Fig. 0, we find that the chemical potential decreases, 
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FIG. 5: (Color online) Circles and squares show the FN-DMC 
energy Edmc for a four-component Fermi gas with N = 16 
and 17 atoms, respectively, interacting through a square well 
potential with range Ro = O.Olaho as a function of the inter- 
species scattering length a s (all interspecies scattering lengths 
are equal). While the N = 16 system has four atoms per 
component, the TV = 17 system has one component with 
five atoms and three components with four atoms. Dotted 
lines show a linear fit to the FN-DMC energies: The slope is 
37.7hco/a ho for N = 16 and A0.7hw/a ho for TV = 17. For com- 
parison, the solid line shows the energy for TV = 16 calculated 
within first order perturbation theory, E — 36hu> + ca s /aho 
with c » 37.1?kj. 



just as the energy, approximately linearly with increas- 
ing \a s \. For scattering lengths more negative than about 
— 0.0450^0 84], the DMC walkers sample — just as in the 
case of the small Bose systems — negative energy config- 
urations. We find that the three-component Fermi gas 
with N = 12 behaves similarly. At the FN-DMC level, 
the instability of multi-component Fermi gases is accom- 
panied by the existence of cluster-like states with neg- 
ative energy. To investigate whether these cluster-like 
states "live" for sufficiently small \a s \, as in the case of 
Bose gases, in a Hilbert space that is approximately or- 
thogonal to the Hilbert space where the gas-like states 
live, we perform a series of VMC calculations. 

Circles and squares in Fig. [5] show the variational en- 
ergy EyMC f° r the three-component Fermi system with 
TV = 12 as a function of b for a s = — 0.05a^ o and —O.laho, 
respectively (as before, b = b a and a s = a a p where 
a ^ 0). The lowest variational energy for a s — — 0.05a^ o 
is obtained for b ~ aho, Evmc — 26.09(4)fiu. The 
FN-DMC energy for this interspecies scattering length, 
Edmc = 26. 08 (4) hoj, agrees with Evmc within error 
bars, indicating that our variational wave function — 
assuming the nodal surface is adequate — captures the 
key physics of the system. For comparison, circles and 
squares in Fig. [7] show the variational energy Evmc for 
the four-component Fermi system with TV = 16 as a func- 
tion of b for a s — — O.Obaho and — Q.Q7aho, respectively. A 
comparison of Figs. [5] and [7] shows that the overall behav- 
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FIG. 6: (Color online) Circles and squares show the varia- 
tional energy Evmc for a three-component Fermi gas with 
TV = 12 atoms interacting through a square well potential 
(range Ro — O.Olaho) with a s = —0.05a,ho and —O.laho (all 
interspecies scattering lengths are equal), respectively, as a 
function of the variational parameter b. Evmc is scaled by 
the energy Emi of the non-interacting system, Eni ~ 27hu). 
Dotted lines connect data points for ease of viewing. 




FIG. 7: (Color online) Circles and squares show the vari- 
ational energy Evmc for a four-component Fermi gas with 
TV = 16 atoms interacting through a square well potential 
(range Ro = 0.01a ho ) with a s = — 0.05a; lo and — 0.07eih o , re- 
spectively, as a function of the variational parameter b. Evmc 
is scaled by the energy Eni of the non-interacting system, 
Emi = 36hu>. Dotted lines connect data points for ease of 
viewing. 



ior of the three- and four-component systems is very sim- 
ilar, and resembles that discussed above for Bose gases: 
For small \a s \, Evmc increases if the Gaussian width b 
is larger than about aho and if b is smaller than about 
aho- For somewhat more negative interspecies scattering 
lengths a s , Ev mc shows a local minimum at b w a^o and 
a local maximum at b w 0.3a^ o , and Evmc becomes neg- 
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FIG. 8: (Color online) Energies for a three-component Fermi 
system as a function of aho/a,s for two resonant interactions. 
Circles and diamonds show the energy of the three-particle 
system, multiplied by a factor of four, interacting through a 
square well potential with range OMa^o and a Gaussian po- 
tential with range 0.007ajj O , respectively; the energies for the 
square well potential are calculated by the FN-DMC approach 
and those for the Gaussian potential by the CG approach. 
Squares and pluses show the FN-DMC energies for TV = 12 
for a square well potential with range O.Ola^o and a Gaussian 
potential with range 0.007ah o , respectively. Solid and dashed 
lines connect data points for ease of viewing. 



ative for yet smaller b. For more negative a s (not shown 
in Figs. [6] and [7|), the energy barrier disappears. We note 
that the variational energies for b values of the order of 
0.3a/j O to 0.6a/j O and a s < — 0.05a/j O show large fluctua- 
tions. For these systems, configurations with rather dif- 
ferent "geometries", and hence with rather different po- 
tential and kinetic energies, are being sampled, leading 
to non-Gaussian distributed energy expectation values. 
This behavior is a consequence of the fact that the varia- 
tional wave function does not provide a good description 
of the true eigenfunction. 

Our MC results for multi-component Fermi gases with 
all resonant interactions and \ > 2 support the physical 
picture developed in Sec. IIIII within the hyperspherical 
framework about how the instability or "implosion" of 
the system arises. Furthermore, the VMC and FN-DMC 
results suggest that the stability and decay of multi- 
component Fermi gases can be attributed to the same 
mechanisms as the stability and decay of Bose gases. At 
the VMC level, the energy barrier arises on length scales 
that are comparable to \a s \ but much larger than the 
range of the two-body potential. This suggests that the 
onset of the instability is determined by the value of the 
scattering length as opposed to the details of the un- 
derlying two-body potential, assuming the range of the 
under-lying two-body potential is sufficiently small (see 
also the discussion in Sec. IIIIBI) . 

We now discuss the behavior of three-component Fermi 
gases with equal masses and equal trapping frequencies in 



which two interspecies scattering lengths, denoted in the 
following by a s , are resonant and the third interspecies 
scattering length is zero. As before, the depths of the 
two-body potentials are adjusted so that the scattering 
lengths take the desired value and so that the two-body 
potentials support no two-body bound state (for a a p = 0, 
the depth is set to zero). 

Motivated by the discussion presented in Sec. IIII Bl for 
the three-particle system we first consider two-body po- 
tentials with very small range. The trimer system in- 
teracting through a square well potential with Rq = 
0.001a^ o supports a state with negative energy for scat- 
tering lengths that are much less negative than the crit- 
ical scattering length of a c « — laho predicted for N = 3 
(see Sec. IIII B|) . Thus, we expect the stability behavior 
of the three-component many-fermion system with two 
resonant interactions interacting through square well po- 
tentials with range Ro = O.OOlaho to be similar to that of 
the all-resonant system discussed above. Indeed, we find 
that the FN-DMC calculations for the N — 12 system ap- 
pear stable for small |a s |, including a region of scattering 
lengths where the three-particle system supports nega- 
tive energy states, but show instabilities for larger \a s \ 
(\a s \ < \a c \). The VMC energy is minimal for b « a/ lo 
and increases for smaller and larger b values. Somewhat 
surprisingly, our VMC calculations do not indicate the 
existence of an energy minimum for small R values when 
a s is close to a c . This is different from the all-resonant 
systems and can possibly be attributed to the variational 
wave functions employed (for very small Rq, a more flex- 
ible functional form may be needed to obtain an energy 
minimum at small R). We find similar results for the 
four-component system with three resonant interactions 
for which dip = a s and a a p — otherwise (we did not 
perform Monte Carlo calculations for the second non- 
trivial configuration with x ~ 1 resonant interactions, for 
which a s = ayi = &23 = and a a p = otherwise). 
Our FN-DMC results for multi-component Fermi gases 
with x ~ 1 resonant interactions (x > 2) and very small 
Ro suggest that these systems become unstable for a cer- 
tain critical negative scattering length. In analogy to the 
Bose system, this seems to suggest that multi-component 
Fermi gases with very small Ro are stable for scattering 
lengths less negative than a c and become mechanically 
unstable for scattering lengths more negative than a c . 

Next, we consider three-component Fermi gases in- 
teracting through two-body potentials with a somewhat 
larger range with a s — a\% — 023 and 031 = 0. Figured] 
shows the energy for the three-particle system, multi- 
plied by a factor of four, as a function of a^ Q /a s for the 
square well potential with Rq = 0.01a^ o (circles) and the 
Gaussian potential with Rq = 0.007a^ o (diamonds). The 
N = 3 energies for the Gaussian potential are calculated 
by the CG approach [52], |H, [H, [H| (this approach is free 
of any assumptions, and its accuracy can be improved 
systematically) and those for the square well potential by 
the FN-DMC approach. For small \a s \, the energies for 
the two different interaction potentials agree quite well 
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TABLE II: Summary of our stability analysis: In the large 
\a a \ limit, trapped x~ com P onen t Fermi gases are predicted 
to be either stable ("S") or unstable ("U"), or to form a 
(possibly stable) gas of trimers ("S(?)GT"). The predictions 
are based on the hyperspherical treatment, which employs 
the bare and the density-dependent pseudo-potential, and a 
many-body Monte Carlo framework. ^Only systems with 
a s — ai2 = ai3 = 014 and a a p — otherwise with very small 
two-body range were investigated. 



X 


hyperspherical 


hyperspherical 


FN-DMC/ 




("bare" PP) 


(den. dep. PP) 


VMC 
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S 


3 (all res.) 


U 
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U 


3 (two res.) 
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U or S(?)GT 


4 (all res.) 
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U 


4 (three res.) 


u 


s 


Tj(i) 



but for large \a s \ deviations are visible. Thus, the sys- 
tem's behavior is, although qualitatively similar, to some 
extent non-universal. The three-body energy becomes 
negative for a s ~ — l.baho and — 3a^ Q for the square well 
and Gaussian potentials, respectively. 

The energies for the three-component Fermi system 
with four atoms per component, i.e., for N = 12, are 
shown in Fig. [8] by squares and pluses for the square 
well and the Gaussian potential, respectively. As in the 
three-particle case, the energies for these two different 
potentials agree well for small \a s \ but behave somewhat 
differently for large \a s \. For both potentials, the a s value 
at which negative energy states appear is similar to that 
for the corresponding three-body system. An analysis 
of the structural properties suggests that the N = 12 
system is made up of weakly-bound three-body clusters. 
Figure [5] shows that the energy for N — 12 is larger 
than four times the energy of the corresponding trimer 
system, suggesting that the many-body system, which 
appears to be made up of weakly-bound trimers, may 
be stable. Our FN-DMC calculations show no evidence 
for the formation of negative energy states consisting of 
clusters that contain more than three atoms. At this 
point it is not clear whether this is a consequence of the 
particular guiding function employed or whether many- 
body clusters consisting of four, five or six atoms with 
negative energy do indeed not exist. 



V. SUMMARY AND CONCLUSION 

This paper investigates the stability of trapped x~ 
component Fermi gases interacting through short-range 
two-body potentials with negative scattering lengths. Ta- 
ble |TT] summarizes the results of our stability analysis. 

One of our main findings is that trapped multi- 
component Fermi gases with more than two components, 
in which all scattering lengths, masses and trapping fre- 



quencies are equal, become unstable for a certain critical 
negative scattering length. This instability, a type of 
Ferminova, is similar to the instability of trapped Bosc 
gases with negative scattering length, sometimes referred 
to as bosenova [25]. The instability predicted here for 
multi-component Fermi systems, \ > 2, emerges within 
two different theoretical frameworks: A hyperspherical 
many-body treatment that employs either the bare or a 
density-dependent zero-range interaction [37|, [HJ , and a 
many-body Monte Carlo treatment, which includes both 
VMC and FN-DMC calculation and employs finite-range 
two-body interactions with a range Rq chosen so that the 
first three-particle state with negative energy appears for 
a s much less negative than a c . 

Both the hyperspherical and the Monte Carlo frame- 
works parameterize the nodal surface of the anti- 
symmetric many-body wave function by the ideal gas 
nodal surface. In the limit of vanishing confinement 
and periodic boundary conditions, this corresponds to a 
wave function that describes a normal state. This paper 
does not investigate wave functions that include pair- 
ing physics or the possibility of phase separation from 
the outset. For three-component systems, for instance, 
a BCS mean-field framework predicts the existence of a 
pairing mechanism in which the atoms of species one and 
two pair (forming a superfluid) while the atoms of species 
three are in a normal state [3(| [H, [H, [35j|. This pair- 
ing mechanism has also been investigated within other 
theoretical frameworks 28, 2j|, [34],[36j]. Although the pa- 
rameter combinations considered in the above cited refer- 
ences is different from those considered in the present pa- 
per, we believe that future Monte Carlo studies of multi- 
component Fermi gases with more than two components 
should investigate the implications of a many-body wave 
function that has pairing physics build in from the outset. 
In particular, it would be interesting to address the ques- 
tion whether a state that contains pairing physics would 
be stable against the collapse investigated in the present 
study. If such a stabilizing mechanism existed this would 
have important consequences for on-going cold atom ex- 
periments. 

To observe the predicted collapse of multi-component 
Fermi gases with all-resonant interactions, we imagine 
the following experiment. First, a stable ^-component 
Fermi gas (x > 2) with vanishing interspecies scattering 
lengths is prepared. We then imagine that all interspecies 
scattering lengths be suddenly tuned to the same large 
negative value. The time scale for this ramp should be 
fast compared to the time scale at which phase sepa- 
ration might occur and also fast compared to the time 
scale at which losses due to spin-flip collisions might oc- 
cur. The attraction between particles would then pro- 
duce an implosion and collapse of the system. Follow- 
ing the implosion, we anticipate that recombination into 
molecules and clusters will immediately follow once the 
system reaches high density, and the resulting energy re- 
lease should detonate a ferminova akin to the bosenova 
that has been observed experimentally and discussed the- 
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oretically [M H H IM Si HI- The rich dynamics 
of soliton formation, akin to that observed for bosonic 
atoms in Ref. [9(| could also conceivably ensue. 

A second main finding of our work concerns three- 
component Fermi systems with two resonant and one 
non-resonant interspecies scattering lengths. The behav- 
ior of these systems depends strongly on the range Rq of 
the underlying two-body potential. For very small i?o, an 
instability similar to that for the three-component sys- 
tem with all-resonant interactions appears to arise, al- 
though at much more negative a s . For larger R Q , the 
many-body system may be made up of weakly-bound 
trimers and may be stable all the way up to unitarity. 
We note that the importance of the range parameter, 
or a three-body parameter, was already pointed out in 
earlier work [3ll . [35l . I3H ] . Future studies should investi- 
gate in more detail whether a Fermi system consisting 
of weakly- or strongly-bound trimers could be realized 
experimentally. 

Finally, we note that multi-component Fermi gases 
have recently also been investigated in one-dimensional 
space [H, |33, E3|- Assuming zero-range interactions, 
many properties of the multi-component system can be 



obtained analytically through the exact Bethe ansatz. 
For negative coupling constant, the ground state of x~ 
component Fermi gases contains clusters that consist 
of x-atoms (one of each species). It is suggested that 
one-dimensional three-component Fermi gases undergo a 
phase transition that is analogous to the quark color su- 
perconductor state [40]. In the future, it would be inter- 
esting to investigate whether such a transition also exists 
in three-dimensional three-component Fermi gases. 

Clearly, more work is needed to better understand the 
behavior of three-dimensional multi-component Fermi 
gases, including the implications of different two-body 
ranges and different nodal surfaces of the guiding func- 
tion employed in the Monte Carlo study. This work in- 
vestigated the behavior of multi-component Fermi gases 
for different scattering lengths and two-body ranges. In 
cases where three-body states with negative energy ex- 
ist, it will be important to determine if the range of the 
two-body parameter is indeed the relavant parameter; it 
may be that the trimer binding energy is a more natural 
quantity. 

We gratefully acknowledge support by the NSF and 
discussions with J. D'Incao. 
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